set.seed(31337)

rm(list = ls())

numiter <- 2500
alpha <- 0.275

setwd("~/Dropbox/complieropt/Polmeth Archive")
source("icsw.R")

library(foreign)
library(rgenoud)
library(gtools)
library(sem)

albertsonlawrence <- read.dta("albertsonlawrence.dta")

attach(albertsonlawrence)

covmat <- cbind(partyid, pnintst, watchnat, readnews, educad, income, gender, white, partyid_d, pnintst_d, watchnat_d, readnews_d, educad_d, income_d, gender_d, race_mv)

partyid[partyid_d == 1] <- mean(partyid[partyid_d == 0])
pnintst[pnintst_d == 1] <- mean(pnintst[pnintst_d == 0])
watchnat[ watchnat_d == 1] <- mean( watchnat[ watchnat_d == 0])
readnews[readnews_d == 1] <- mean(readnews[readnews_d == 0])
educad[educad_d == 1] <- mean(educad[educad_d == 0])
gender[gender_d == 1] <- mean(gender[gender_d == 0])
income[income_d == 1] <- mean(income[income_d == 0])
white[race_mv == 1] <- mean(white[race_mv == 0])

covmatH <- cbind(partyid,pnintst, watchnat,educad,readnews,gender,income,white)

Y = infopro2 # alternative DV: support3
X = watchpro
Z = conditn

cscoreout <- complier.score(X,Z,covmatH,genoud=FALSE)

summary(cscoreout$C.score)

cscore <- cscoreout$C.score

N <- length(Y)

qcscore <- quantile(cscore, 1/N^alpha)
cscore[cscore < qcscore] <- qcscore

Ymis1 <- is.na(infopro2)

IPWweight1<- 1/(1-predict(glm(Ymis1~Z+covmatH,family=binomial(link="probit")),type="response"))
IPWweight1[Ymis1] <- 0

outputTSLS <- tsls.wfit(cbind(1,watchpro,covmatH),infopro2,cbind(1,conditn,covmatH),w=IPWweight1)

outputICSW <- tsls.wfit(cbind(1,watchpro,covmatH),infopro2,cbind(1,conditn,covmatH),w=IPWweight1/cscore)

outputTSLS$coefficients
outputICSW$coefficients

####

Ymis2 <- is.na(support3)

IPWweight2<- 1/(1-predict(glm(Ymis2~Z+covmatH,family=binomial(link="probit")),type="response"))
IPWweight2[Ymis2] <- 0

outputTSLS.2 <- tsls.wfit(cbind(1,watchpro,covmatH),support3,cbind(1,conditn,covmatH),w=IPWweight2)

outputICSW.2 <- tsls.wfit(cbind(1,watchpro,covmatH),support3,cbind(1,conditn,covmatH),w=IPWweight2/cscore)

outputTSLS.2$coefficients
outputICSW.2$coefficients

ests <- matrix(NA,numiter,40)

###

for (i in 1:numiter) {

N <- length(watchpro)

sampleN <- c(sample(c(1:N),replace=TRUE))

#[conditn==1],replace=TRUE),sample(c(1:N)[conditn==0],replace=TRUE))

IPWweight1BS<- 1/(1-predict(glm(Ymis1[sampleN]~Z[sampleN]+covmatH[sampleN,],family=binomial(link="probit"),maxit=2000),type="response"))
IPWweight1BS[Ymis1[sampleN]] <- 0

IPWweight2BS<- 1/(1-predict(glm(Ymis2[sampleN]~Z[sampleN]+covmatH[sampleN,],family=binomial(link="probit"),maxit=2000),type="response"))
IPWweight2BS[Ymis2[sampleN]] <- 0

cscoreoutBS <- complier.score(watchpro[sampleN],conditn[sampleN],covmatH[sampleN,],genoud=FALSE)

cscoreBS <- cscoreoutBS$C.score
qcscoreBS <- quantile(cscoreBS, 1/N^alpha)
cscoreBS[cscoreBS < qcscoreBS] <- qcscoreBS

outputTSLSbs <- tsls.wfit(cbind(1,watchpro,covmatH)[sampleN,],infopro2[sampleN],cbind(1,conditn,covmatH)[sampleN,],w=IPWweight1BS)

outputICSWbs <- tsls.wfit(cbind(1,watchpro,covmatH)[sampleN,],infopro2[sampleN],cbind(1,conditn,covmatH)[sampleN,],w=IPWweight1BS/cscoreBS)

outputTSLSbs.2 <- tsls.wfit(cbind(1,watchpro,covmatH)[sampleN,],support3[sampleN],cbind(1,conditn,covmatH)[sampleN,],w=IPWweight2BS)

outputICSWbs.2 <- tsls.wfit(cbind(1,watchpro,covmatH)[sampleN,],support3[sampleN],cbind(1,conditn,covmatH)[sampleN,],w=IPWweight2BS/cscoreBS)


ests[i,] <- c(outputTSLSbs$coefficients,
outputICSWbs$coefficients,outputTSLSbs.2$coefficients,
outputICSWbs.2$coefficients)

}

colMeans(ests)
apply(ests,2,sd)

mean(ests[,12]-ests[,2])
quantile(ests[,12]-ests[,2],c(0.025,0.975))
mean(ests[,12]-ests[,2] <= 0)

mean(ests[,32]-ests[,22])
quantile(ests[,32]-ests[,22],c(0.025,0.975))
mean(ests[,32]-ests[,22] <= 0)

######

cor(cbind(cscore,covmatH))
cor(cbind(cscoreout$A.score,covmatH))

###

perm.test <- function(D,Z,Xmat,cscoreout,numperm=50,genoud=FALSE) {
	test.dist <- rep(NA,numperm)
	
	for (i in 1:numperm) {
		Xmats <- Xmat[sample(1:length(D)),]
		cout <- complier.score(D,Z,Xmats,genoud=genoud)
		test.dist[i] <- sum((D - cout$C.score*Z - cout$A.score)^2)
	}
	
	test.stat <- sum((D - cscoreout$C.score*Z - cscoreout$A.score)^2)
	p.value <- mean(test.dist <= test.stat)
	output <- list(test.dist,test.stat,p.value)
	names(output) <- c("test.dist","test.stat","p.value")
	return(output)
}


permout <- perm.test(X,Z,covmatH,cscoreout,numperm=numiter)

permout$p.value

plot(density(permout$test.dist))
lines(c(permout$test.stat,permout$test.stat),c(-100,100))

save.image(file="albertsonlawrenceresults.RData")

detach(albertsonlawrence)

load("albertsonlawrenceresults.RData")

#### cscore matrix

xtable(round(cor(cbind(cscore,covmatH)),digits=2))

#### cscore p-value

permout$p.value

#### Results: Information

# 2SLS Estimate
outputTSLS$coefficients

# ICSW Estimate
outputICSW$coefficients

# Bootstrap p-value Difference in Coefficients
quantile(ests[,12]-ests[,2],probs=c(0.025,0.05,0.95,0.975))

#### Results: Support

# 2SLS Estimate
outputTSLS.2$coefficients

# ICSW Estimate
outputICSW.2$coefficients

# Bootstrap p-value Difference in Coefficients
quantile(ests[,32]-ests[,22],probs=c(0.025,0.05,0.95,0.975))

round(cbind(outputTSLS$coefficients,outputICSW$coefficients,outputTSLS.2$coefficients,outputICSW.2$coefficients),digits=2)

mat1 <- round(cbind(apply(ests,2,quantile,c(0.025))[1:10],apply(ests,2,quantile,c(0.025))[11:20],apply(ests,2,quantile,c(0.025))[21:30],apply(ests,2,quantile,c(0.025))[31:40]),digits=2)

mat2 <- round(cbind(apply(ests,2,quantile,c(0.975))[1:10],apply(ests,2,quantile,c(0.975))[11:20],apply(ests,2,quantile,c(0.975))[21:30],apply(ests,2,quantile,c(0.975))[310]),digits=2)

summary(ests)


quartz(width=4,height=4)

par(mfrow = c(1, 1), pty = "s",family="Gill Sans MT",font.main=1,cex=.7)

plot(density(permout$test.dist),xlim=c(min(permout$test.dist),max(permout$test.dist)+ 1),main="Albertson and Lawrence (2009) Permutation Test",xlab="SSR")
lines(c(permout$test.stat,permout$test.stat),c(-100,100))
lines(c(-100,10000),c(0,0))
